!###############################################
 module energy
!###############################################
   implicit none
   integer :: ndim_T, m_T, nnz_T
   real (kind=8) :: ntime,timeref
   integer, dimension (:), allocatable :: rowind_T,columns_T
   integer, dimension (:), allocatable :: irn_T, jcn_T
   real (kind=8), dimension (:), allocatable :: rhs_T,sol_T,nnzval_T
   real (kind=8), dimension (:,:), allocatable :: mtrx_T,al_T,b_T
   real (kind=8), dimension(:), allocatable :: T_z,v_z,ux_z,uz_z,eta_z,Phi_z,Gam_z,RaH_z, &
                                               alpha_z,kap_z,depl_z,Tsol_z

!*******************
  contains
!*******************

!============================================================
 subroutine addiff
! Solution of the advection-diffusion equation of temperature
! and/or advection of composition
!============================================================
   use share; use domain; use markeraux; use flow; use linsolver; use iniconds
   implicit none
   integer :: i,j,ir,rlen
   real (kind=8) :: a,dt_fix,dt_Cnr
   character (len=6) :: cha
!============================================================
  
!-------------------------------
! Generate 2D grid and markers
!-------------------------------
   call gridgen

!--------------------------------------------
! Read input variables from yacc.inp
!--------------------------------------------
   call read_input_file

!---------------------------------------------------------
! Allocate the variables for solving the Stokes equation
!---------------------------------------------------------
   call stokes_alloc

!-----------------------------------------------
! Setup boundary conditions for Stokes equation
!-----------------------------------------------
   call set_bc_flow

!--------------------------------------------------------
! Boundary conditions for the energy equation
! bcb_T = 0 bottom temperature is fixed at Tbott
!         1 bottom heat flux is fixed at qbott
! at lateral boundaries zero heat flux is prescribed
! at the top fixed temperature Ttop is prescribed
!--------------------------------------------------------
   if(bcB_T.eq.0) then
     Ttop=0d0
     Tbott=1d0     
   else
     Ttop=0d0
     qbott=0d0
   endif  

!-------------------------------------------------------
! Memory allocations for advection-diffusion equation
!-------------------------------------------------------
   call addiff_alloc

!----------------------------------------
! Set internal heating RaH as the array
!----------------------------------------
   RaH=RaHval

!-----------------------------------------
! Set time-step index and initial time
!-----------------------------------------
   if(restart.eq.0) then
     timin=0
     ntime=0d0
     timeref=0d0
   else
     open(1,file=trim(outdir)//'grd_restart.dat')
     read(1,*) cha,cha,timin,cha,ntime
     close(1)
     timax=timin+timax
     timeref=ntime
   endif

!========================  Henri's code to output velocity ===================================
!  if(binout.eq.1) then
!    rlen=8*((nx+2)*(nz+2)+1)    !Record length for nodal points & time in single precision    
!    OPEN(101,FILE='vx_bin.dat',FORM='UNFORMATTED',STATUS='REPLACE',ACCESS='DIRECT',RECL=rlen)
!    OPEN(102,FILE='vz_bin.dat',FORM='UNFORMATTED',STATUS='REPLACE',ACCESS='DIRECT',RECL=rlen)
!  endif
!=============================================================================================

   ir=0
!*********************
time: do ti=timin,timax
!*********************

!============================================
! index for record length of binary output
   ir=ir+1 
!============================================

!+++++++++++++++++++++
   if(ti.eq.timin) then
!+++++++++++++++++++++
   
!---------------------
! Initial conditions
!---------------------

! Read initial conditions from file...
   if(restart.eq.1) then
     
     print *, 'Restart from previous run'

     open(1,file=trim(outdir)//'grd_restart.dat')
     read(1,*) cha,cha,timin,cha,ntime
     do i=0,nx+1
       do j=0,nz+1
         read(1,*) a,a,a,a,a,Tnew(i,j),eta(i,j),a,Gam(i,j),RaH(i,j)
       end do
     end do
     if(bcB_T.eq.0) Tbott=Tnew(1,0)
     close(1)

     if(markers.eq.1) then

       do i=1,Nmk
         if(comp_mk.eq.1) then
           Gmk(i)=0d0
           call grid2mk(xmk(i),zmk(i),Gam,a)
           Gmk(i)=a
         endif

         if(heats_mk.eq.1) then
           RaHmk(i)=0d0
           call grid2mk(xmk(i),zmk(i),RaH,a)
           RaHmk(i)=a
         endif

         if(visc_mk.eq.1) then
           etamk(i)=0d0
           call grid2mk(xmk(i),zmk(i),eta,a)
           etamk(i)=a
         endif

       end do
     endif

   elseif(restart.eq.0) then

! ...otherwise define the initial distribution of
! temperature and/or...
     if(therm.eq.1) then
       if(debug.eq.1) print *, 'Initialize temperature'
       call initial_temp
       if(debug.eq.1) print *, '   Done'
     endif

!...other fields defined on markers
     if(markers.eq.1) then
       if(debug.eq.1) print *, 'Initialize distribution of markers properties'
       call initial_markers
       if(debug.eq.1) print *, '   Done'
     endif

!...plus melt-related stuff
     if(melt.eq.1) then
       if(debug.eq.1) print *, 'Initialize depletion, solidus and liquidus'
       call initial_melt
       if(debug.eq.1) print *, '   Done'
     endif

   endif

!+++++++++
   else
!+++++++++
   
! Thermal convection only
   if(therm.eq.1.and.markers.eq.0) then
     if(debug.eq.1) print *, 'Solve energy'
     call T_SLCN_OS
     if(debug.eq.1) print *, '   Done'

! Compositional convection only
   elseif(therm.eq.0.and.markers.gt.0) then
     if(debug.eq.1) print *, 'Markers advection and interpolation'
     call mkadvect
     if(debug.eq.1) print *, '   Done'
   
! Thermo-chemical convection
   elseif(therm.eq.1.and.markers.gt.0) then
     if(debug.eq.1) print *, 'Solve energy'
     call T_SLCN_OS
     if(debug.eq.1) print *, '   Done'

! Melting
     if(melt.eq.1) then
       if(debug.eq.1) print *, 'Markers advection and interpolation'
       call calc_melt
       if(debug.eq.1) print *, '   Done'
     endif

     if(debug.eq.1) print *, 'Markers advection and interpolation'
     call mkadvect
     if(debug.eq.1) print *, '   Done'

   endif

! Decay of heat sources
   if(decayQ.eq.1.and.heats_mk.eq.0) then
     if(debug.eq.1) print *, 'Decay of heat sources on grid'
     call heat_decay
     if(debug.eq.1) print *, '   Done'
   endif

! Calculate new CMB temperature 
   if(core_cooling.eq.1) then
     if(debug.eq.1) print *, 'Core cooling'
     call update_Tbott
     if(debug.eq.1) print *, '   Done'
   endif

!+++++++++++
   endif
!+++++++++++

!-------------------------------------------------------------
! Solve the Stokes equation with the new temperature 
! and concentration and choose the time-step dt for a given
! Courant number 'Cnr'
!-------------------------------------------------------------
   uold=u
   vold=v
   if(debug.eq.1) print *, 'Solve momentum'
   call stokes
   if(debug.eq.1) print *, '   Done'

!--------------------------------------------------
! Time-extrapolated velocity used to integrate the 
! trajectory equation backwards in time
!--------------------------------------------------
   if(ti.eq.0.or.ext_fac.eq.1d0) then 
     uextr=u
     vextr=v
   else
     uextr = ext_fac*u + (1d0-ext_fac)*uold
     vextr = ext_fac*v + (1d0-ext_fac)*vold
   endif

!--------------------------------------------
! Time stepping according to Courant number
!--------------------------------------------
   if(Cnr.ne.0d0) then
     dt=Cnr*min(minval(delx)/maxval(dabs(u)),minval(delz)/maxval(dabs(v)))
     dt_Cnr=dt
   endif

   if((init_dt.gt.0d0).and.(ti.eq.0).and.(Cnr.ne.0d0).and.(restart.eq.0)) then
     dt=init_dt
     dt_fix=2d0*init_dt
   elseif((init_dt.gt.0d0).and.(ti.gt.0).and.(dt_fix.le.dt_Cnr).and.(restart.eq.0)) then
     dt=dt_fix
     dt_fix=2d0*dt_fix
   endif

! In case of restart, the solution is recomputed and the time is not advanced
   if(restart.eq.1.and.ti.eq.timin) cycle

   ntime=ntime+dt

!------------------------------------------------
! Compute averages, profiles and store output
!------------------------------------------------
   if(debug.eq.1) print *, 'Calculate averages'
   call timeseries(ti)
   if(debug.eq.1) print *, '   Done'

   if((ti.eq.timin)                                                  &
       .or.(stp_ascii.gt.0.and.mod(ti,stp_ascii).eq.0)               &
       .or.(stp_ascii.eq.0.and.ntime-timeref.ge.time_ascii)) then 

     if(debug.eq.1) print *, 'Calculate vertical profiles'
     call calc_profiles
     if(debug.eq.1) print *, '   Done'
     
     if(debug.eq.1) print *, 'Write output files'
     call write_output
     if(debug.eq.1) print *, '   Done'

   endif       

   if(ntime.gt.time_max) stop

!**************** 
   end do time
!**************** 


 end subroutine addiff
 

!=============================================
  subroutine addiff_alloc
! Various allocations needed in addiff
!=============================================
   use share; use domain
   allocate(Told(0:nx+1,0:nz+1),Tnew(0:nx+1,0:nz+1),drho(0:nx+1,0:nz+1),              & 
            T_z(0:nz+1),v_z(0:nz+1),ux_z(0:nz+1),uz_z(0:nz+1),eta_z(0:nz+1),          &
            Phi_z(0:nz+1),Gam_z(0:nz+1),RaH_z(0:nz+1),alpha_z(0:nz+1),kap_z(0:nz+1),  &
            Tsol_z(0:nz+1),depl_z(0:nz+1))
   Told=0d0; Tnew=0d0 
   allocate(Tpr(0:nx+1,0:nz+1))     
   Tpr=0d0

!---------------------------------------------------------------------------------------
! ndim_T => integer.  Size of the matrix for implicit solution of the heat equation
! nnz_T => integer. Number of non-zero elements of mtrx_T
! nnzval_T => integer, array(nnz). Non-zero elements stored row after row
! columns_T => integer, array(nnz). Column index of nonzero elements
! rowind_T => integer, array(ndim+1). Array containing the incremental numbers 
!                                   of nonzero elements that appear in a new row
!---------------------------------------------------------------------------------------
   if(bcB_T.eq.0) then
     ndim_T=(nx+2)*nz
   else  
     ndim_T=(nx+2)*(nz+1)
   endif

   allocate(rhs_T(ndim_T),b_T(ndim_T,1),sol_T(ndim_T))

   nnz_T=5*(nx+2)*(nz+2)
   allocate(nnzval_T(nnz_T),columns_T(nnz_T),rowind_T(ndim_T+1))
   allocate(irn_T(nnz_T),jcn_T(nnz_T))

   allocate(Gam(0:nx+1,0:nz+1),Gam_old(0:nx+1,0:nz+1),eta(0:nx+1,0:nz+1),eta_old(0:nx+1,0:nz+1),   &
            RaH(0:nx+1,0:nz+1),RaH_old(0:nx+1,0:nz+1),depl(0:nx+1,0:nz+1),depl_old(0:nx+1,0:nz+1), &
            Tsol0(0:nz+1),Tsol(0:nx+1,0:nz+1),Tsol_old(0:nx+1,0:nz+1),Tliq(0:nz+1))
            
   if(markers.gt.0) then    
     if(comp_mk.eq.1) allocate(Gmk(Nmk),dGmk(Nmk))
     if(heats_mk.eq.1) allocate(RaHmk(Nmk))
     if(visc_mk.eq.1) allocate(etamk(Nmk))
     if(depl_mk.eq.1) allocate(deplmk(Nmk))
     allocate(Tmk(Nmk),umk(Nmk),vmk(Nmk))
   endif

 end subroutine addiff_alloc

!========================================
 subroutine mkadvect
!========================================
   use share; use domain; use markeraux; use flow 
   implicit none
   integer :: i,j,l,rko
   integer, dimension(4) :: ai,aj
   real (kind=8) :: xm,zm,dm,xmnew,zmnew
   real (kind=8), dimension(4) :: w
   real (kind=8), dimension(0:nx+1,0:nz+1) :: pom
!==============================================================
 
!--------------------------------------------------
! Passive advection of markers via Runge-Kutta
!--------------------------------------------------
   rko=4
   do i=1,Nmk
     xm=xmk(i)
     zm=zmk(i)
     select case(rko)
     case(2) !2nd order Runge-Kutta
     call rk2(xm,zm,xmnew,zmnew)
     case(4) !4th order Runge-Kutta
     call rk4(xm,zm,xmnew,zmnew)
     end select
     xmk(i)=xmnew
     zmk(i)=zmnew
   end do     

!---------------------------------------------------------------------------------------------------------
! Change markers properties if necessary
! Only implemented for magma ocean overturn with density differences associated with phase transitions
!---------------------------------------------------------------------------------------------------------
   if(icd.eq.5) then

     do i=1,Nmk
       dm=dimz-zmk(i)

       if((mktype(i).eq.2).and.(dm.le.0.578d0)) then
         mktype(i)=1
         Gmk(i)=Gmk(i)-dGmk(i)
       endif

       if((mktype(i).eq.1).and.(dm.ge.0.578d0)) then
         mktype(i)=2
         Gmk(i)=Gmk(i)+dGmk(i)
       endif

     end do

   endif

!------------------------------------------------------------------------------
! Calculate new fields by interpolating the marker properties
! from the advected markers onto the grid
!------------------------------------------------------------------------------
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   if(comp_mk.eq.1.or.heats_mk.eq.1.or.visc_mk.eq.1.or.depl_mk.eq.1) then
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

   if(comp_mk.eq.1) Gam=0d0
   if(heats_mk.eq.1) RaH=0d0
   if(visc_mk.eq.1) eta=0d0
   if(depl_mk.eq.1) depl=0d0

! Decay of heat sources on markers
   if(decayQ.eq.1.and.heats_mk.eq.1) then
     if(debug.eq.1) print *, 'Decay of heat sources on markers'
     do l=1,Nmk
       RaHmk(l)=RaHmk(l)*exp(-decay_lambda*dt)
     end do
     if(debug.eq.1) print *, 'Done'
   endif

   pom=0d0
   do l=1,Nmk
     xm=xmk(l)
     zm=zmk(l)

!  Compute the weights of all markers found in the 4 cells
!  that surround each grid point 
    call Wmarker4(xm,zm,ai,aj,w)
    do i=1,4 
      pom(ai(i),aj(i))=pom(ai(i),aj(i))+w(i)
      if(comp_mk.eq.1) Gam(ai(i),aj(i))=Gam(ai(i),aj(i))+Gmk(l)*w(i)
      if(heats_mk.eq.1) RaH(ai(i),aj(i))=RaH(ai(i),aj(i))+RaHmk(l)*w(i)
      if(visc_mk.eq.1) eta(ai(i),aj(i))=eta(ai(i),aj(i))+etamk(l)*w(i)
      if(depl_mk.eq.1) depl(ai(i),aj(i))=depl(ai(i),aj(i))+deplmk(l)*w(i)
    end do  

   end do

   if(comp_mk.eq.1) then
     Gam=Gam/pom

     do i=0,nx+1
       do j=0,nz+1
         if (pom(i,j).eq.0d0) then
            if(debug.eq.1) write(*,*) "Grid node ",i,j,"has zero weights. Use C-value from previous time step."
            Gam(i,j)=Gam_old(i,j)
         endif
       end do
     end do
     Gam_old=Gam     

   endif

   if(heats_mk.eq.1) then
     RaH=RaH/pom

     do i=0,nx+1
       do j=0,nz+1
         if (pom(i,j).eq.0d0) then
            if(debug.eq.1) write(*,*) "Grid node ",i,j,"has zero weights. Use RaH-value from previous time step."
            RaH(i,j)=RaH_old(i,j)
         endif
       end do
     end do
     RaH_old=RaH     

   endif

   if(visc_mk.eq.1) then
     eta=eta/pom

     do i=0,nx+1
       do j=0,nz+1
         if (pom(i,j).eq.0d0) then
            if(debug.eq.1) write(*,*) "Grid node ",i,j,"has zero weights. Use RaH-value from previous time step."
            eta(i,j)=eta_old(i,j)
         endif
       end do
     end do
     eta_old=eta     

   endif

   if(depl_mk.eq.1) then
     depl=depl/pom

     do i=0,nx+1
       do j=0,nz+1
         Tsol(i,j)=Tsol(i,j) + (depl(i,j)-depl_old(i,j))*(Tliq(j)-Tsol0(j))
       end do
     end do
 
     do i=0,nx+1
       do j=0,nz+1
         if (pom(i,j).eq.0d0) then
            if(debug.eq.1) write(*,*) "Grid node ",i,j,"has zero weights. Use depl-value from previous time step."
            depl(i,j)=depl_old(i,j)
         endif
       end do
     end do

   endif

!++++++++++++
   endif
!++++++++++++

 end subroutine mkadvect

 
!======================================================================
 subroutine T_SLCN_OS
! Semi Lagrangian Crank Nicolson Operator Splitting 
! solution to the advection-diffusion equation
! Spiegelman & Katz, G^3 (2006)
!
! First the advection equation is solved, then diffusion and (where
! applies) reaction, are applied to the advected temperature field
!======================================================================
   use share; use domain; use markeraux; use flow
   implicit none
   integer :: i,j,l,k,i1,j1,interp,jmin,rko,ndtau
   real (kind=8) :: xst,zst,xold,zold,dtau,T,xl,xu,zl,zu,ddum
   real (kind=8) :: z,rho_r,drho_r,T_r,z0,Cl,Tph,fph,lath
   real (kind=8), dimension(:), allocatable :: ff   
   real (kind=8), dimension(4) :: y,y1,y2,y12
   real (kind=8), dimension(0:nx+1,0:nz+1) :: y2a
   real (kind=8), parameter :: ee=1d-12
   logical :: rk,substep
   fph(z)=0.5d0*(1d0+dtanh((z-z0-Cl*(T-Tph))/dph))
!======================================================================

   Told=Tnew
!---------------------------------------------------------------------
! According to the extrapolated velocity field (uextr,vextr), 
! find the take-off point, interpolate there Told and find the 
! solution Tpr of the advection equation only
! xst: x of the take-off point
! zst: z of the take-off point
! interp: order of the interpolation
!       =1 for bilinear interpolation
!       =2 for bicubic interpolation
!       =3 for bicubic spline interpolation
! rk: how to find the take-off point
!     =.true. via 2nd or 4th order Runge-Kutta
!     =.false. via single-step Euler method
! rko: order of the Runge-Kutta integrator
! substep: if .true. a time substepping is used to find the take-off
!          point via Runge-Kutta 
! dtau: time substep for trajectory equation
! ndtau: number of substeps to be taken within dt
!---------------------------------------------------------------------
   interp=2
   rk=.true.
   rko=4
   substep=.true.
   if(Cnr.gt.1d0) then
     ndtau=int(2d0*Cnr)
     dtau=dt/ndtau
   else
     substep=.false.
   endif

   if(bcB_T.eq.0) then
     jmin=1
   else
     jmin=0
   endif

! Construct second derivatives y2a for spline interpolation
   if(interp.eq.3) call splie2(xgrid,zgrid,Told,nx+1,nz+1,y2a)

   do i=0,nx+1
     do j=jmin,nz     
     
!---------------------------------------
!    Find the take-off point (xst,zst)
!---------------------------------------
        if(rk) then

          if(rko.eq.2) then

            if(substep) then      

              xold=xgrid(i); zold=zgrid(j)
              do k=1,ndtau       
                call rk2bt_substep(xold,zold,dtau,xst,zst)
                xold=xst; zold=zst
              end do

            else
              call rk2bt(i,j,xst,zst)
            endif

          elseif(rko.eq.4) then

            if(substep) then 

              xold=xgrid(i); zold=zgrid(j)
              do k=1,ndtau       
                call rk4bt_substep(xold,zold,dtau,xst,zst)
                xold=xst; zold=zst
              end do

            else
              call rk4bt(i,j,xst,zst)
            endif

          endif

        else

           xst=xgrid(i)-uextr(i,j)*dt
           zst=zgrid(j)-vextr(i,j)*dt

         endif

         if(xst.lt.0d0.or.xst.gt.dimx+ee.or.  &
            zst.lt.0d0.or.zst.gt.dimz+ee) then
            print *, 'take off point out of domain'
            print *, 'x(i)=',xgrid(i),'x*=',xst,'z*=',zst,'z(j)=',zgrid(j)
         endif

       select case(interp)

         case(1)
!-----------------------------------
!       Bilinear interpolation         	 
!-----------------------------------
           call grid2mk(xst,zst,Told,T)

         case(2)
!-----------------------------------
!       Bicubic interpolation...         	 
!-----------------------------------
           i1=bisec(xgrid,xst)
           j1=bisec(zgrid,zst)
           call derivs(i1,j1,Told,y,y1,y2,y12)
           xl=xgrid(i1); xu=xgrid(i1+1)
           zl=zgrid(j1); zu=zgrid(j1+1)
           call bcuint(y,y1,y2,y12,xl,xu,zl,zu,xst,zst,T,ddum,ddum)

         case(3)
!-----------------------------------
!       Spline interpolation         	 
!-----------------------------------
           call splin2(xgrid,zgrid,Told,y2a,nx+1,nz+1,xst,zst,T)

       end select
       Tpr(i,j)=T

     end do
   end do

   if(bcB_T.eq.0) then
     Tpr(:,0)=Tbott
     Tpr(:,nz+1)=Ttop
   else
     Tpr(:,nz+1)=Ttop
   endif

!-----------------------------------------------------------------------------
! For all grid points, apply to 'Tpr' the operator (I + dt/2(nabla^2 - Di v)) 
! and setup the rhs to be used as source for the diffusion.
! I: identity 
! dt: time step 
! nabla^2: discrete Laplace operator
! Di: dissipation number
! v: vertical flow
!-----------------------------------------------------------------------------
   if(nphases.ne.0) allocate(ff(nphases))

   l=0
   do i=0,nx+1
     do j=jmin,nz
       l=l+1      

       z=dimz-zgrid(j)
       T=Tpr(i,j)

       if(cpr.eq.1) then
         call refstate(zgrid(j),rho_r,drho_r,T_r)
       else
         rho_r=1d0
         drho_r=1d0
         T_r=0d0
       endif

! Latent heat 
       lath=0d0
       if(nphases.ne.0) then
         
         do k=1,nphases
           z0=phase_z(k)
           Cl=phase_Cl(k)
           Tph=phase_T(k)
           ff(k)=fph(z)

           lath=lath+dt/2/rho_r**2*phase_Rb(k)/Ra*Di*v(i,j)*phase_Cl(k)*2d0/dph*(ff(k)-ff(k)**2)*Tpr(i,j)   &
                    +dt/rho_r**2*phase_Rb(k)/Ra*Di*v(i,j)*phase_Cl(k)*2d0/dph*(ff(k)-ff(k)**2)*Tsurf   
         end do
!         lath=0d0

       endif

       if(j.eq.0) then
!     Bottom boundary nodes: prescribed heat flux 'qbott'       

         if(i.eq.0) then
           rhs_T(l)=Tpr(i,j)+dt/2d0/rho_r*                                                                  &
                    (4d0/(delx(i+1)+delx(i+1))*((kap(i,j)+kap(i+1,j))/2d0*(Tpr(i+1,j)-Tpr(i,j))/delx(i+1))  &
                    +4/(delz(j+1)+delz(j+1))*(kap(i,j)+kap(i,j+1))/2d0*(Tpr(i,j+1)-Tpr(i,j))/delz(j+1)      &
                    -4/(delz(j+1)+delz(j+1))*qbott)                                                         &
                    -dt/2*Di*alpha(i,j)*v(i,j)*Tpr(i,j)-dt*Di*alpha(i,j)*v(i,j)*Tsurf                       &
                    +dt/rho_r*Di/Ra*Phi(i,j)+dt*RaH(i,j)/Ra+lath


         elseif(i.eq.nx+1) then
           rhs_T(l)=Tpr(i,j)+dt/2/rho_r*                                                                    &
                    (-4/(delx(i)+delx(i))*((kap(i,j)+kap(i-1,j))/2d0*(Tpr(i,j)-Tpr(i-1,j))/delx(i))         &
                    +4/(delz(j+1)+delz(j+1))*(kap(i,j)+kap(i,j+1))/2d0*(Tpr(i,j+1)-Tpr(i,j))/delz(j+1)      &
                    -4/(delz(j+1)+delz(j+1))*qbott)                                                         &
                    -dt/2*Di*alpha(i,j)*v(i,j)*Tpr(i,j)-dt*Di*alpha(i,j)*v(i,j)*Tsurf                       &
                    +dt/rho_r*Di/Ra*Phi(i,j)+dt*RaH(i,j)/Ra+lath

         else
           rhs_T(l)=Tpr(i,j)+dt/2/rho_r*                                                                    &
                    (2/(delx(i)+delx(i+1))*((kap(i,j)+kap(i+1,j))/2d0*(Tpr(i+1,j)-Tpr(i,j))/delx(i+1)       &
                                           -(kap(i,j)+kap(i-1,j))/2d0*(Tpr(i,j)-Tpr(i-1,j))/delx(i))        &
                    +4/(delz(j+1)+delz(j+1))*(kap(i,j)+kap(i,j+1))/2d0*(Tpr(i,j+1)-Tpr(i,j))/delz(j+1)      &
                    -4/(delz(j+1)+delz(j+1))*qbott)                                                         &
                    -dt/2*Di*alpha(i,j)*v(i,j)*Tpr(i,j)-dt*Di*alpha(i,j)*v(i,j)*Tsurf                       &
                    +dt/rho_r*Di/Ra*Phi(i,j)+dt*RaH(i,j)/Ra+lath

         endif

       else
         
         if(i.eq.0) then
!     Left boundary nodes: zero heat flux
           rhs_T(l)=Tpr(i,j)+dt/2d0/rho_r*                                                                  &
                    (4d0/(delx(i+1)+delx(i+1))*((kap(i,j)+kap(i+1,j))/2d0*(Tpr(i+1,j)-Tpr(i,j))/delx(i+1))  &
                    +2d0/(delz(j)+delz(j+1))*((kap(i,j)+kap(i,j+1))/2d0*(Tpr(i,j+1)-Tpr(i,j))/delz(j+1)     &
                                             -(kap(i,j)+kap(i,j-1))/2d0*(Tpr(i,j)-Tpr(i,j-1))/delz(j)))     &
                    -dt/2*Di*alpha(i,j)*v(i,j)*Tpr(i,j)-dt*Di*alpha(i,j)*v(i,j)*Tsurf                       &
                    +dt/rho_r*Di/Ra*Phi(i,j)+dt*RaH(i,j)/Ra+lath

         elseif(i.eq.nx+1) then
!     Right boundary nodes: zero heat flux
           rhs_T(l)=Tpr(i,j)+dt/2d0/rho_r*                                                                  &
                    (-4d0/(delx(i)+delx(i))*((kap(i,j)+kap(i-1,j))/2d0*(Tpr(i,j)-Tpr(i-1,j))/delx(i))       &
                    +2d0/(delz(j)+delz(j+1))*((kap(i,j)+kap(i,j+1))/2d0*(Tpr(i,j+1)-Tpr(i,j))/delz(j+1)     & 
                                           -(kap(i,j)+kap(i,j-1))/2d0*(Tpr(i,j)-Tpr(i,j-1))/delz(j)))       &
                    -dt/2*Di*alpha(i,j)*v(i,j)*Tpr(i,j)-dt*Di*alpha(i,j)*v(i,j)*Tsurf                       &
                    +dt/rho_r*Di/Ra*Phi(i,j)+dt*RaH(i,j)/Ra+lath

         else  
!     Internal nodes       
           rhs_T(l)=Tpr(i,j)+dt/2/rho_r*                                                                    &
                    (2/(delx(i)+delx(i+1))*((kap(i,j)+kap(i+1,j))/2d0*(Tpr(i+1,j)-Tpr(i,j))/delx(i+1)       &
                                   -(kap(i,j)+kap(i-1,j))/2d0*(Tpr(i,j)-Tpr(i-1,j))/delx(i))                &
                    +2/(delz(j)+delz(j+1))*((kap(i,j)+kap(i,j+1))/2d0*(Tpr(i,j+1)-Tpr(i,j))/delz(j+1)       &
                                           -(kap(i,j)+kap(i,j-1))/2d0*(Tpr(i,j)-Tpr(i,j-1))/delz(j)))       &
                    -dt/2*Di*alpha(i,j)*v(i,j)*Tpr(i,j)-dt*Di*alpha(i,j)*v(i,j)*Tsurf                       &
                    +dt/rho_r*Di/Ra*Phi(i,j)+dt*RaH(i,j)/Ra+lath

         endif

       endif 

     end do
   end do  

   if(nphases.ne.0) deallocate(ff)

   call diffusion_SLCN

 end subroutine T_SLCN_OS


!==========================================================
 subroutine diffusion_SLCN
! Solution of the linear system of the diffusion equation 
! via Crank-Nicolson scheme
!==========================================================
   use share; use domain; use linsolver; use flow
   implicit none
   integer :: i,j,l,k,jmin,innz
   real (kind=8) :: z,T,rho_r,drho_r,T_r,fph,lath,z0,Cl,Tph
   real (kind=8), dimension(:), allocatable :: ff
   fph(z)=0.5d0*(1d0+dtanh((z-z0-Cl*(T-Tph))/dph))
   real (kind=8), dimension(:,:), allocatable :: e_matrix
   character (len=32) :: varfmt
!======================================================================

   if(nphases.ne.0) allocate(ff(nphases))

!-----------------------------------------------------------------
! Setup compactly the matrix resulting from the operator 
! (I - dt/2(nabla^2 - Di v)) applied to the unknown temperature
!-----------------------------------------------------------------
   sol_T=0d0
   nnzval_T=0d0
   columns_T=0
   rowind_T=0
   innz=0

   if(bcB_T.eq.0) then
     jmin=1
   else
     jmin=0
   endif

   l=0
   do i=0,nx+1
     do j=jmin,nz
       l=l+1

       z=dimz-zgrid(j)
       T=Told(i,j)

       if(cpr.eq.1) then
         z=zgrid(j)
         call refstate(z,rho_r,drho_r,T_r)
       else
         rho_r=1d0
         drho_r=1d0
       endif

! Latent heat 
       if(nphases.ne.0) then
         
         lath=0d0
         do k=1,nphases
           z0=phase_z(k)
           Cl=phase_Cl(k)
           Tph=phase_T(k)
           ff(k)=fph(z)

           lath=lath+dt/2/rho_r**2*phase_Rb(k)/Ra*Di*v(i,j)*phase_Cl(k)*2d0/dph*(ff(k)-ff(k)**2)   
         end do

       endif
       
!    ***********************
       if(bcB_T.eq.0) then ! prescribed bottom temperature
!    ***********************

       rowind_T(l)=innz+1   

       if(i.eq.0) then 

         if(j.ne.1) then
         innz=innz+1
         nnzval_T(innz)=-(dt/2)/rho_r*2/(delz(j)+delz(j+1))*(kap(i,j)+kap(i,j-1))/2d0*1/delz(j)    !T(i,j-1)
         columns_T(innz)=l-1
         endif

         innz=innz+1
         nnzval_T(innz)=1d0+(dt/2)/rho_r*(4/(delx(i+1)+delx(i+1))*(kap(i,j)+kap(i+1,j))/2d0*(1/delx(i+1))    & !T(i,j)
                                 +2/(delz(j)+delz(j+1))*((kap(i,j)+kap(i,j+1))/2d0*1/delz(j+1)               & 
                                                        +(kap(i,j)+kap(i,j-1))/2d0*1/delz(j)))               &  
                                 +(dt/2)*Di*alpha(i,j)*v(i,j)                                                &
                                 +lath
         columns_T(innz)=l   

         if(j.ne.nz) then
         innz=innz+1
         nnzval_T(innz)=-(dt/2)/rho_r*2/(delz(j)+delz(j+1))*(kap(i,j)+kap(i,j+1))/2d0*1/delz(j+1)   !T(i,j+1)
         columns_T(innz)=l+1
         endif

         innz=innz+1
         nnzval_T(innz)=-(dt/2)/rho_r*4/(delx(i+1)+delx(i+1))*(kap(i,j)+kap(i+1,j))/2d0*1/delx(i+1) !T(i+1,j)
         columns_T(innz)=l+nz

       elseif(i.eq.nx+1) then 

         innz=innz+1
         nnzval_T(innz)=-(dt/2)/rho_r*4/(delx(i)+delx(i))*(kap(i-1,j)+kap(i,j))/2d0*1/delx(i)   !T(i-1,j)
         columns_T(innz)=l-nz

         if(j.ne.1) then
         innz=innz+1
         nnzval_T(innz)=-(dt/2)/rho_r*2/(delz(j)+delz(j+1))*(kap(i,j-1)+kap(i,j))/2d0*1/delz(j)    !T(i,j-1)
         columns_T(innz)=l-1
         endif  

         innz=innz+1
         nnzval_T(innz)=1+(dt/2)/rho_r*(4/(delx(i)+delx(i))*(kap(i-1,j)+kap(i,j))/2d0*(1/delx(i))       & !T(i,j)
                              +2/(delz(j)+delz(j+1))*((kap(i,j+1)+kap(i,j))/2d0*1/delz(j+1)+            &
                                                      (kap(i,j-1)+kap(i,j))/2d0*1/delz(j)))             & 
                              +(dt/2)*Di*alpha(i,j)*v(i,j)                                              &
                              +lath
         columns_T(innz)=l

         if(j.ne.nz) then
         innz=innz+1
         nnzval_T(innz)=-(dt/2)/rho_r*2/(delz(j)+delz(j+1))*(kap(i,j)+kap(i,j+1))/2d0*1/delz(j+1)  !T(i,j+1)
         columns_T(innz)=l+1
         endif

       else

         innz=innz+1
         nnzval_T(innz)=-(dt/2)/rho_r*2/(delx(i)+delx(i+1))*(kap(i-1,j)+kap(i,j))/2d0*1/delx(i)   !T(i-1,j)
         columns_T(innz)=l-nz

         if(j.ne.1) then
         innz=innz+1
         nnzval_T(innz)=-(dt/2)/rho_r*2/(delz(j)+delz(j+1))*(kap(i,j-1)+kap(i,j))/2d0*1/delz(j)    !T(i,j-1)
         columns_T(innz)=l-1
         endif

         innz=innz+1
         nnzval_T(innz)=1+(dt/2)/rho_r*(2/(delx(i)+delx(i+1))*((kap(i+1,j)+kap(i,j))/2d0*1/delx(i+1)  &  !T(i,j)
                                                              +(kap(i-1,j)+kap(i,j))/2d0*1/delx(i))   &
                                       +2/(delz(j)+delz(j+1))*((kap(i,j+1)+kap(i,j))/2d0*1/delz(j+1)  &
                                                              +(kap(i,j-1)+kap(i,j))/2d0*1/delz(j)))  & 
                                       +(dt/2)*Di*alpha(i,j)*v(i,j)                                   &
                                       +lath
         columns_T(innz)=l

         if(j.ne.nz) then
         innz=innz+1
         nnzval_T(innz)=-(dt/2)/rho_r*2/(delz(j)+delz(j+1))*(kap(i,j+1)+kap(i,j))/2d0*1/delz(j+1)  !T(i,j+1)
         columns_T(innz)=l+1
         endif

         innz=innz+1
         nnzval_T(innz)=-(dt/2)/rho_r*2/(delx(i)+delx(i+1))*(kap(i+1,j)+kap(i,j))/2d0*1/delx(i+1) !T(i+1,j)
         columns_T(innz)=l+nz
      
       endif

!    ***************************
       elseif(bcB_T.eq.1) then ! prescribed bottom heat flux
!    ***************************

       rowind_T(l)=innz+1   

       if(j.eq.0) then

         if(i.eq.0) then 
           innz=innz+1
           nnzval_T(innz)=1+(dt/2)/rho_r*(4/(delx(i+1)+delx(i+1))*(kap(i+1,j)+kap(i,j))/2d0*1/delx(i+1)     & !T(i,j)
                                +4/(delz(j+1)+delz(j+1))*(kap(i,j+1)+kap(i,j))/2d0*1/delz(j+1))             & 
                                +(dt/2)*Di*alpha(i,j)*v(i,j)                                                &   
                                +lath
           columns_T(innz)=l

           innz=innz+1
           nnzval_T(innz)=-(dt/2)/rho_r*4/(delz(j+1)+delz(j+1))*(kap(i,j+1)+kap(i,j))/2d0*1/delz(j+1)      !T(i,j+1)
           columns_T(innz)=l+1

           innz=innz+1
           nnzval_T(innz)=-(dt/2)/rho_r*4/(delx(i+1)+delx(i+1))*(kap(i+1,j)+kap(i,j))/2d0*1/delx(i+1)         !T(i+1,j)
           columns_T(innz)=l+(nz+1)

         elseif(i.eq.nx+1) then

           innz=innz+1
           nnzval_T(innz)=-(dt/2)/rho_r*4/(delx(i)+delx(i))*(kap(i-1,j)+kap(i,j))/2d0*1/delx(i)               !T(i-1,j)
           columns_T(innz)=l-(nz+1)

           innz=innz+1
           nnzval_T(innz)=1+(dt/2)/rho_r*(4/(delx(i)+delx(i))*(kap(i-1,j)+kap(i,j))/2d0*1/delx(i)    & !T(i,j)
                                +4/(delz(j+1)+delz(j+1))*(kap(i,j+1)+kap(i,j))/2d0*1/delz(j+1))      & 
                                +(dt/2)*Di*alpha(i,j)*v(i,j)                                         & 
                                +lath
           columns_T(innz)=l

           innz=innz+1
           nnzval_T(innz)=-(dt/2)/rho_r*4/(delz(j+1)+delz(j+1))*(kap(i,j+1)+kap(i,j))/2d0*1/delz(j+1)      !T(i,j+1)
           columns_T(innz)=l+1

         else
         
           innz=innz+1
           nnzval_T(innz)=-(dt/2)/rho_r*2/(delx(i)+delx(i+1))*(kap(i-1,j)+kap(i,j))/2d0*1/delx(i)             !T(i-1,j)
           columns_T(innz)=l-(nz+1)

           innz=innz+1
           nnzval_T(innz)=1+(dt/2)/rho_r*(2/(delx(i)+delx(i+1))*((kap(i+1,j)+kap(i,j))/2d0*1/delx(i+1)   & !T(i,j)
                                                                +(kap(i-1,j)+kap(i,j))/2d0*1/delx(i))    & 
                                         +4/(delz(j+1)+delz(j+1))*(kap(i,j+1)+kap(i,j))/2d0*1/delz(j+1)) & 
                                         +(dt/2)*Di*alpha(i,j)*v(i,j)                                    &
                                         +lath
           columns_T(innz)=l

           innz=innz+1
           nnzval_T(innz)=-(dt/2)/rho_r*4/(delz(j+1)+delz(j+1))*(kap(i,j+1)+kap(i,j))/2d0*1/delz(j+1)    !T(i,j+1)
           columns_T(innz)=l+1

           innz=innz+1
           nnzval_T(innz)=-(dt/2)/rho_r*2/(delx(i)+delx(i+1))*(kap(i+1,j)+kap(i,j))/2d0*1/delx(i+1)      !T(i+1,j)	 
           columns_T(innz)=l+(nz+1)
         endif

       else

         if(i.eq.0) then 

           innz=innz+1
           nnzval_T(innz)=-(dt/2)/rho_r*2/(delz(j)+delz(j+1))*(kap(i,j-1)+kap(i,j))/2d0*1/delz(j)     !T(i,j-1)
           columns_T(innz)=l-1

           innz=innz+1
           nnzval_T(innz)=1+(dt/2)/rho_r*(4/(delx(i+1)+delx(i+1))*(kap(i+1,j)+kap(i,j))/2d0*(1/delx(i+1))    &   !T(i,j)
                                +2/(delz(j)+delz(j+1))*((kap(i,j+1)+kap(i,j))/2d0*1/delz(j+1)                &
                                                       +(kap(i,j-1)+kap(i,j))/2d0*1/delz(j)))                & 
                                +(dt/2)*Di*alpha(i,j)*v(i,j)                                                &
                                +lath
           columns_T(innz)=l

           if(j.ne.nz) then
           innz=innz+1
           nnzval_T(innz)=-(dt/2)/rho_r*2/(delz(j)+delz(j+1))*(kap(i,j+1)+kap(i,j))/2d0*1/delz(j+1)    !T(i,j+1)
           columns_T(innz)=l+1
           endif

           innz=innz+1
           nnzval_T(innz)=-(dt/2)/rho_r*4/(delx(i+1)+delx(i+1))*(kap(i+1,j)+kap(i,j))/2d0*1/delx(i+1)   !T(i+1,j)
           columns_T(innz)=l+(nz+1)

         elseif(i.eq.nx+1) then 

           innz=innz+1
           nnzval_T(innz)=-(dt/2)/rho_r*4/(delx(i)+delx(i))*(kap(i-1,j)+kap(i,j))/2d0*1/delx(i)         !T(i-1,j)
           columns_T(innz)=l-(nz+1)

           innz=innz+1
           nnzval_T(innz)=-(dt/2)/rho_r*2/(delz(j)+delz(j+1))*(kap(i,j-1)+kap(i,j))/2d0*1/delz(j)       !T(i,j-1)
           columns_T(innz)=l-1

           innz=innz+1
           nnzval_T(innz)=1+(dt/2)/rho_r*(4/(delx(i)+delx(i))*(kap(i-1,j)+kap(i,j))/2d0*(1/delx(i))     &       !T(i,j)
                                +2/(delz(j)+delz(j+1))*((kap(i,j+1)+kap(i,j))/2d0*1/delz(j+1)           &
                                                       +(kap(i,j-1)+kap(i,j))/2d0*1/delz(j)))           & 
                                +(dt/2)*Di*alpha(i,j)*v(i,j)                                            &
                                +lath
           columns_T(innz)=l

           if(j.ne.nz) then
           innz=innz+1
           nnzval_T(innz)=-(dt/2)/rho_r*2/(delz(j)+delz(j+1))*(kap(i,j+1)+kap(i,j))/2d0*1/delz(j+1)    !T(i,j+1)
           columns_T(innz)=l+1
           endif

         else

           innz=innz+1
           nnzval_T(innz)=-(dt/2)/rho_r*2/(delx(i)+delx(i+1))*(kap(i-1,j)+kap(i,j))/2d0*1/delx(i)       !T(i-1,j)
           columns_T(innz)=l-(nz+1)

           innz=innz+1
           nnzval_T(innz)=-(dt/2)/rho_r*2/(delz(j)+delz(j+1))*(kap(i,j-1)+kap(i,j))/2d0*1/delz(j)       !T(i,j-1)
           columns_T(innz)=l-1

           innz=innz+1
           nnzval_T(innz)=1+(dt/2)/rho_r*(2/(delx(i)+delx(i+1))*((kap(i+1,j)+kap(i,j))/2d0*1/delx(i+1)  &   !T(i,j)
                                                                +(kap(i-1,j)+kap(i,j))/2d0*1/delx(i))   &
                                         +2/(delz(j)+delz(j+1))*((kap(i,j+1)+kap(i,j))/2d0*1/delz(j+1)  &
                                                                +(kap(i,j-1)+kap(i,j))/2d0*1/delz(j)))  & 
                                         +(dt/2)*Di*alpha(i,j)*v(i,j)                                           &
                                         +lath
           columns_T(innz)=l

           if(j.ne.nz) then
           innz=innz+1
           nnzval_T(innz)=-(dt/2)/rho_r*2/(delz(j)+delz(j+1))*(kap(i,j+1)+kap(i,j))/2d0*1/delz(j+1)    !T(i,j+1)
           columns_T(innz)=l+1
           endif

           innz=innz+1
           nnzval_T(innz)=-(dt/2)/rho_r*2/(delx(i)+delx(i+1))*(kap(i+1,j)+kap(i,j))/2d0*1/delx(i+1)    !T(i+1,j)
           columns_T(innz)=l+(nz+1)

         endif
       endif    

!    ********
       endif
!    ********

!--------------------------------
! Setup the right hand side
!--------------------------------

!    ***********************
       if(bcB_T.eq.0) then 
!    ***********************

         if(j.eq.1) then 
!    Bottom temperature fixed at Tbott
           rhs_T(l)=rhs_T(l)+(dt/2)/rho_r*2/(delz(j)+delz(j+1))/delz(j)*(kap(i,j-1)+kap(i,j))/2d0*Tbott
         elseif(j.eq.nz) then 
!    Surface temperature fixed at Ttop
           rhs_T(l)=rhs_T(l)+(dt/2)/rho_r*2/(delz(j)+delz(j+1))/delz(j+1)*(kap(i,j+1)+kap(i,j))/2d0*Ttop
         endif

!    ***************************
       elseif(bcB_T.eq.1) then 
!    ***************************

         if(j.eq.0) then 
!    Bottom heat flux fixed at qbott
           rhs_T(l)=rhs_T(l)-(dt/2)/rho_r*4/(delz(j+1)+delz(j+1))*qbott
         elseif(j.eq.nz) then 
!    Surface temperature fixed at Ttop
           rhs_T(l)=rhs_T(l)+(dt/2)/rho_r*2/(delz(j)+delz(j+1))/delz(j+1)*(kap(i,j+1)+kap(i,j))/2d0*Ttop
         endif

!    **********
       endif
!    **********

     end do 
   end do
   
   if(nphases.ne.0) deallocate(ff)
    
   rowind_T(ndim_T+1)=innz+1

!--------------------------------
! Solution of the linear system
!--------------------------------
   sol_T=0d0
   b_T(:,1)=rhs_T

   call mtrx_idxs(ndim_T,columns_T,rowind_T,irn_T,jcn_T)

! Retrieve and write on file the whole energy matrix (only used for debug/development) 
   if(write_e_matrix.eq.1) then
     print *, 'Store energy matrix in e_matrix.dat'      

     allocate(e_matrix(ndim_T,ndim_T))      
     call mtrx_restore(nnzval_T,columns_T,rowind_T,e_matrix,ndim_T)

     open(1,file='e_matrix.dat')
     write(varfmt,99) ndim_T
     do i=1,ndim_T
        write(1,varfmt) (e_matrix(i,j), j=1,ndim_T) 
     end do
     deallocate(e_matrix)
99   format ('(',I0,'E12.3)')
   endif  

!----------------------------------------
!     Solution with MUMPS
!----------------------------------------
   mumps_par%JOB = -1
   mumps_par%SYM = 0
   mumps_par%PAR = 1
   call dmumps(mumps_par)
   mumps_par%ICNTL(1) = 0
   mumps_par%ICNTL(2) = 0
   mumps_par%ICNTL(3) = minfo
   mumps_par%ICNTL(4) = 2
   mumps_par%ICNTL(7) = mord
   mumps_par%ICNTL(23) = 4000
   mumps_par%N = ndim_T
   mumps_par%NZ = nnz_T      
   allocate(mumps_par%IRN ( mumps_par%NZ ) )
   allocate(mumps_par%JCN ( mumps_par%NZ ) )
   allocate(mumps_par%A( mumps_par%NZ ) )
   allocate(mumps_par%RHS ( mumps_par%N  ) )

   do i = 1, mumps_par%NZ
     mumps_par%IRN(i) = irn_T(i)
     mumps_par%JCN(i) = jcn_T(i)
     mumps_par%A(i) = nnzval_T(i)
   end do

   do i = 1, mumps_par%N
     mumps_par%RHS(i) = rhs_T(i)
   end do

!  Call package for solution
   mumps_par%JOB = 6
   CALL DMUMPS(mumps_par)

!  Store the mumps solution 
   sol_T = mumps_par%RHS

!  Deallocate user data
   deallocate( mumps_par%IRN )
   deallocate( mumps_par%JCN )
   deallocate( mumps_par%A   )
   deallocate( mumps_par%RHS )

!  Destroy the instance (deallocate internal data structures)
   mumps_par%JOB = -2
   CALL DMUMPS(mumps_par)

!---------------------------------------------------
! Store the solution into the new temperature Tnew
!---------------------------------------------------
   l=0
   do i=0,nx+1
     do j=jmin,nz       
       l=l+1
       Tnew(i,j)=sol_T(l)
     end do
   end do

!----------------------------------------------------
! Set boundary temperatures
!----------------------------------------------------

   if(bcB_T.eq.0) then 
     Tnew(:,0)=Tbott
     Tnew(:,nz+1)=Ttop
   elseif(bcB_T.eq.1) then 
     Tnew(:,nz+1)=Ttop
   endif
   
 end subroutine diffusion_SLCN

!=================================================================
 subroutine heat_decay
! Decay of heat sources on the grid accoring to a single 
! average decay constant
!=================================================================
   use share; use domain
   integer :: i,j
!=================================================================

   do i=0,nx+1
     do j=0,nz+1
       RaH(i,j)=RaH(i,j)*exp(-decay_lambda*dt)
     end do
   end do

 end subroutine heat_decay

!===========================================================================
 subroutine update_Tbott
! Core cooling: update the bottom temperature according to the equation:
! dT(z=0)/dt = -3*rho_m*cp_m/(rho_c*cp_c*Rad_c)* dT(z=0)/dz
!===========================================================================
   use share; use domain
   integer :: i
   real (kind=8) :: Tz0,Tz1
!=================================================================
   
! Compute the the temperature gradient at the bottom boundary
   Tz0=0d0; Tz1=0d0
   do i=1,nx+1
     Tz0=Tz0 + delx(i)/2d0*(Tnew(i,0)+Tnew(i-1,0))
     Tz1=Tz1 + delx(i)/2d0*(Tnew(i,1)+Tnew(i-1,1))
   end do
   Tz0=Tz0/dimx
   Tz1=Tz1/dimx

   Tbott = Tbott + dt*3d0*cooling_factor*(Tz1-Tz0)/delz(1)

 end subroutine update_Tbott

!=========================================================================
 subroutine calc_melt
! For all grid-points where T>Tsol, update the temperature, calculate
! the depletion and set the solidus back to the new temperature
!=========================================================================
   use share; use domain; use markeraux; use flow
   implicit none
   integer :: i,j
   real (kind=8) :: val
!======================================================================
   
   do i=0,nx+1
     do j=0,nz+1
       
       if(Tnew(i,j).gt.Tsol(i,j).and.Tnew(i,j).lt.Tliq(j)) then

!      Correct the temperature
         Tnew(i,j) = (Tnew(i,j) + deltaS*Tsol(i,j)*(Tsurf+Tsol(i,j))/(Tliq(j)-Tsol0(j))) &
                    /(1d0 + deltaS*(Tsurf+Tsol(i,j))/(Tliq(j)-Tsol0(j)))

!      Calculate the depletion       
         depl(i,j) = depl(i,j) + (Tnew(i,j)-Tsol(i,j))/(Tliq(j)-Tsol0(j))

!      Update the solidus
         Tsol(i,j) = Tnew(i,j) 
 
       endif

     end do
   end do
   depl_old=depl

! Interpolate the new depletion back on markers	 
   do i=1,Nmk
     deplmk(i)=0d0
     call grid2mk(xmk(i),zmk(i),depl,val)
     deplmk(i)=val
   end do

 end subroutine calc_melt


!=================================================================
 subroutine derivs(k,l,f,y,y1,y2,y12)   
! Given a function tabulated at grid values 'f', and given
! a location (k,l) on the grid, returns the value 'y', the
! x-derivative 'y1', the z-derivative 'y2' and the mixed second
! derivative 'y12' of the tabulated function at 4 cell corners 
! counted couterclockwise starting from (k,l) 
!=================================================================
   use share; use domain
   implicit none
   integer, intent(in) :: k,l
   real (kind=8), dimension(:), intent(out) :: y,y1,y2,y12
   real (kind=8), dimension(0:,0:), intent(in) :: f
!=================================================================
   
   y(1)=f(k,l)
   y(2)=f(k+1,l)
   y(3)=f(k+1,l+1)
   y(4)=f(k,l+1)

! First derivative along x
   if(k.eq.0) then
     y1(1)=(f(1,l)-f(0,l))/delx(1)
     y1(2)=(f(2,l)-f(0,l))/(delx(2)+delx(1))
     y1(3)=(f(2,l+1)-f(0,l+1))/(delx(2)+delx(1))
     y1(4)=(f(1,l+1)-f(0,l+1))/delx(1)
   elseif(k.eq.nx) then
     y1(1)=(f(nx+1,l)-f(nx-1,l))/(delx(nx+1)+delx(nx))
     y1(2)=(f(nx+1,l)-f(nx,l))/delx(nx+1)
     y1(3)=(f(nx+1,l+1)-f(nx,l+1))/delx(nx+1)
     y1(4)=(f(nx+1,l+1)-f(nx-1,l+1))/(delx(nx+1)+delx(nx))
   else  
     y1(1)=(f(k+1,l)-f(k-1,l))/(xgrid(k+1)-xgrid(k-1))
     y1(2)=(f(k+2,l)-f(k,l))/(xgrid(k+2)-xgrid(k))
     y1(3)=(f(k+2,l+1)-f(k,l+1))/(xgrid(k+2)-xgrid(k))
     y1(4)=(f(k+1,l+1)-f(k-1,l+1))/(xgrid(k+1)-xgrid(k-1))
   endif
  
! First derivative along z 
   if(l.eq.0) then
     y2(1)=(f(k,l+1)-f(k,l))/delz(1)
     y2(2)=(f(k+1,l+1)-f(k+1,l))/delz(1)
     y2(3)=(f(k+1,l+2)-f(k+1,l))/(zgrid(l+2)-zgrid(l))
     y2(4)=(f(k,l+2)-f(k,l))/(zgrid(l+2)-zgrid(l))
   elseif(l.eq.nz) then
     y2(1)=(f(k,l+1)-f(k,l-1))/(zgrid(l+1)-zgrid(l-1))
     y2(2)=(f(k+1,l+1)-f(k+1,l-1))/(zgrid(l+1)-zgrid(l-1))
     y2(3)=(f(k+1,l+1)-f(k+1,l))/delz(nz+1)
     y2(4)=(f(k,l+1)-f(k,l))/delz(nz+1)
   else
     y2(1)=(f(k,l+1)-f(k,l-1))/(zgrid(l+1)-zgrid(l-1))
     y2(2)=(f(k+1,l+1)-f(k+1,l-1))/(zgrid(l+1)-zgrid(l-1))
     y2(3)=(f(k+1,l+2)-f(k+1,l))/(zgrid(l+2)-zgrid(l))
     y2(4)=(f(k,l+2)-f(k,l))/(zgrid(l+2)-zgrid(l))
   endif

! Mixed second derivatives
   if(k.eq.0.and.l.eq.0) then
     y12(1)=(f(1,1)-f(1,0)-f(0,1)+f(0,0))  &
            /delx(1)/delz(1)
     y12(2)=(f(2,1)-f(2,0)-f(0,1)+f(0,0))   &
            /(xgrid(2)-xgrid(0))/delz(1)
     y12(3)=(f(2,2)-f(2,0)-f(0,2)+f(0,0)) &
            /(xgrid(2)-xgrid(0))/(zgrid(2)-zgrid(0))
     y12(4)=(f(1,2)-f(1,0)-f(0,2)+f(0,0))  &
            /delx(1)/(delz(1)+delz(2))
   elseif(k.eq.nx.and.l.eq.0) then
     y12(1)=(f(k+1,l+1)-f(k+1,l)-f(k-1,l+1)+f(k-1,l)) &
          /(xgrid(k+1)-xgrid(k-1))/delz(1)
     y12(2)=(f(nx+1,1)-f(nx+1,0)-f(nx,1)+f(nx,0)) &
            /delx(nx+1)/delz(1)
     y12(3)=(f(nx+1,2)-f(nx+1,0)-f(nx,2)+f(nx,0)) &
            /delx(nx+1)/(delz(1)+delz(2))
     y12(4)=(f(k+1,l+2)-f(k+1,l)-f(k-1,l+2)+f(k-1,l)) &
            /(xgrid(k+1)-xgrid(k-1))/(zgrid(l+2)-zgrid(l))
   elseif(k.eq.0.and.l.eq.nz) then    
     y12(1)=(f(1,nz+1)-f(1,nz-1)-f(0,nz+1)+f(0,nz-1)) &
            /delx(1)/(delz(nz)+delz(nz+1))
     y12(2)=(f(k+2,l+1)-f(k+2,l-1)-f(k,l+1)+f(k,l-1)) &
            /(xgrid(k+2)-xgrid(k))/(zgrid(l+1)-zgrid(l-1))
     y12(3)=(f(k+2,l+1)-f(k+2,l)-f(k,l+1)+f(k,l)) &
            /(xgrid(k+2)-xgrid(k))/delz(nz+1)
     y12(4)=(f(1,nz+1)-f(1,nz)-f(0,nz+1)+f(0,nz)) &
            /delx(1)/delz(nz+1)
   elseif(k.eq.nx.and.l.eq.nz) then   
     y12(1)=(f(k+1,l+1)-f(k+1,l-1)-f(k-1,l+1)+f(k-1,l-1)) &
            /(xgrid(k+1)-xgrid(k-1))/(zgrid(l+1)-zgrid(l-1))
     y12(2)=(f(nx+1,nz+1)-f(nx+1,nz-1)-f(nx,nz+1)+f(nx,nz-1)) &
            /delx(nx+1)/(delz(nz)+delz(nz+1))
     y12(3)=(f(nx+1,nz+1)-f(nx+1,nz)-f(nx,nz+1)+f(nx,nz)) &
            /delx(nx+1)/delz(nz+1)
     y12(4)=(f(k+1,l+1)-f(k+1,l)-f(k-1,l+1)+f(k-1,l)) &
            /(xgrid(k+1)-xgrid(k-1))/delz(nz+1)
   elseif(k.eq.0) then
     y12(1)=(f(1,l+1)-f(1,l-1)-f(0,l+1)+f(0,l-1)) &
            /delx(1)/(delz(l)+delz(l+1))
     y12(2)=(f(k+2,l+1)-f(k+2,l-1)-f(k,l+1)+f(k,l-1)) &
            /(xgrid(k+2)-xgrid(k))/(zgrid(l+1)-zgrid(l-1))
     y12(3)=(f(k+2,l+2)-f(k+2,l)-f(k,l+2)+f(k,l)) &
            /(xgrid(k+2)-xgrid(k))/(zgrid(l+2)-zgrid(l))
     y12(4)=(f(1,l+2)-f(1,l)-f(0,l+2)+f(0,l)) &
            /delx(1)/(delz(l+1)+delz(l+2))
   elseif(k.eq.nx) then
     y12(1)=(f(k+1,l+1)-f(k+1,l-1)-f(k-1,l+1)+f(k-1,l-1)) &
            /(xgrid(k+1)-xgrid(k-1))/(zgrid(l+1)-zgrid(l-1))
     y12(2)=(f(nx+1,l+1)-f(nx+1,l-1)-f(nx,l+1)+f(nx,l-1)) &
            /delx(nx+1)/(delz(l+1)+delz(l))
     y12(3)=(f(nx+1,l+2)-f(nx+1,l)-f(nx,l+2)+f(nx,l)) &
            /delx(nx+1)/(delz(l+2)+delz(l+1))
     y12(4)=(f(k+1,l+2)-f(k+1,l)-f(k-1,l+2)+f(k-1,l)) &
            /(xgrid(k+1)-xgrid(k-1))/(zgrid(l+2)-zgrid(l))
   elseif(l.eq.0) then
     y12(1)=(f(k+1,l+1)-f(k+1,l)-f(k-1,l+1)+f(k-1,l)) &
            /(xgrid(k+1)-xgrid(k-1))/delz(1)
     y12(2)=(f(k+2,l+1)-f(k+2,l)-f(k,l+1)+f(k,l)) &
            /(xgrid(k+2)-xgrid(k))/delz(1)
     y12(3)=(f(k+2,l+2)-f(k+2,l)-f(k,l+2)+f(k,l)) &
            /(xgrid(k+2)-xgrid(k))/(zgrid(l+2)-zgrid(l))
     y12(4)=(f(k+1,l+2)-f(k+1,l)-f(k-1,l+2)+f(k-1,l)) &
            /(xgrid(k+1)-xgrid(k-1))/(zgrid(l+2)-zgrid(l))
   elseif(l.eq.nz) then
     y12(1)=(f(k+1,l+1)-f(k+1,l-1)-f(k-1,l+1)+f(k-1,l-1)) &
            /(xgrid(k+1)-xgrid(k-1))/(zgrid(l+1)-zgrid(l-1))
     y12(2)=(f(k+2,l+1)-f(k+2,l-1)-f(k,l+1)+f(k,l-1)) &
            /(xgrid(k+2)-xgrid(k))/(zgrid(l+1)-zgrid(l-1))
     y12(3)=(f(k+2,l+1)-f(k+2,l)-f(k,l+1)+f(k,l)) &
            /(xgrid(k+2)-xgrid(k))/delz(nz+1)
     y12(4)=(f(k+1,l+1)-f(k+1,l)-f(k-1,l+1)+f(k-1,l)) &
            /(xgrid(k+1)-xgrid(k-1))/delz(nz+1)
   else
     y12(1)=(f(k+1,l+1)-f(k+1,l-1)-f(k-1,l+1)+f(k-1,l-1)) &
            /(xgrid(k+1)-xgrid(k-1))/(zgrid(l+1)-zgrid(l-1))
     y12(2)=(f(k+2,l+1)-f(k+2,l-1)-f(k,l+1)+f(k,l-1)) &
            /(xgrid(k+2)-xgrid(k))/(zgrid(l+1)-zgrid(l-1))
     y12(3)=(f(k+2,l+2)-f(k+2,l)-f(k,l+2)+f(k,l)) &
            /(xgrid(k+2)-xgrid(k))/(zgrid(l+2)-zgrid(l))
     y12(4)=(f(k+1,l+2)-f(k+1,l)-f(k-1,l+2)+f(k-1,l)) &
            /(xgrid(k+1)-xgrid(k-1))/(zgrid(l+2)-zgrid(l))
   endif

 end subroutine derivs 


!=================================================
      subroutine splie2(x1a,x2a,ya,m,n,y2a)
!=================================================
      integer :: m,n,j,k
      integer, parameter :: NN=2000
      real (kind=8) :: x1a(0:m),x2a(0:n),y2a(0:m,0:n),ya(0:m,0:n)
      real (kind=8) :: y2tmp(0:NN),ytmp(0:NN)

      do  j=0,m
        do  k=0,n
          ytmp(k)=ya(j,k)
        end do
        call spline(x2a,ytmp,n,1d30,1d30,y2tmp)
        do k=0,n
          y2a(j,k)=y2tmp(k)
        end do
      end do

      end subroutine splie2

!=================================================
      subroutine spline(x,y,n,yp1,ypn,y2)
!=================================================
      integer :: n,i,k
      integer, parameter :: nmax=2000
      real (kind=8) :: yp1,ypn,x(0:n),y(0:n),y2(0:n)
      real (kind=8) :: p,qn,sig,un,u(0:nmax)

      if (yp1.gt..99d30) then
        y2(1)=0.
        u(1)=0.
      else
        y2(1)=-0.5
        u(1)=(3./(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
      endif
      do i=1,n-1
        sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
        p=sig*y2(i-1)+2.
        y2(i)=(sig-1.)/p
        u(i)=(6.*((y(i+1)-y(i))/(x(i+1)-x(i))-(y(i)-y(i-1))/(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*u(i-1))/p
      end do
      if (ypn.gt..99d30) then
        qn=0.
        un=0.
      else
        qn=0.5
        un=(3./(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
      endif
      y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.)
      do k=n-1,0,-1
        y2(k)=y2(k)*y2(k+1)+u(k)
      end do

      end subroutine spline

!=================================================
      subroutine splin2(x1a,x2a,ya,y2a,m,n,x1,x2,y)
!=================================================
      integer :: m,n,j,k
      integer, parameter  :: nn=2000
      real (kind=8) :: x1,x2,y,x1a(0:m),x2a(0:n),y2a(0:m,0:n),ya(0:m,0:n)
      real (kind=8) :: y2tmp(0:nn),ytmp(0:nn),yytmp(0:nn)

      do j=0,m
        do k=0,n
          ytmp(k)=ya(j,k)
          y2tmp(k)=y2a(j,k)
        end do
        call splint(x2a,ytmp,y2tmp,n,x2,yytmp(j))
      end do
      call spline(x1a,yytmp,m,1d30,1d30,y2tmp)
      call splint(x1a,yytmp,y2tmp,m,x1,y)

      end subroutine splin2

!=================================================
      subroutine splint(xa,ya,y2a,n,x,y)
!=================================================
      use share
      integer :: n,khi,klo
      real (kind=8) :: x,y,xa(0:n),y2a(0:n),ya(0:n)
      real (kind=8) :: a,b,h

      klo=bisec(xa,x)
      khi=klo+1
      h=xa(khi)-xa(klo)
      if (h.eq.0d0) print *, 'bad xa input in splint'
      a=(xa(khi)-x)/h
      b=(x-xa(klo))/h
      y=a*ya(klo)+b*ya(khi)+((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))*(h**2)/6.

      end subroutine splint


!=================================================
        subroutine bcucof(y,y1,y2,y12,d1,d2,c)
!=================================================
        implicit none
        real(kind=8), intent(in) :: d1,d2
        real(kind=8), dimension(4), intent(in) :: y,y1,y2,y12
        real(kind=8), dimension(4,4), intent(out) :: c
        real(kind=8), dimension(16) :: x
        real(kind=8), dimension(16,16) :: wt
        data wt /1,0,-3,2,4*0,-3,0,9,-6,2,0,-6,4,&
                8*0,3,0,-9,6,-2,0,6,-4,10*0,9,-6,2*0,-6,4,2*0,3,-2,6*0,-9,6,&
                2*0,6,-4,4*0,1,0,-3,2,-2,0,6,-4,1,0,-3,2,8*0,-1,0,3,-2,1,0,-3,&
                2,10*0,-3,2,2*0,3,-2,6*0,3,-2,2*0,-6,4,2*0,3,-2,0,1,-2,1,5*0,&
                -3,6,-3,0,2,-4,2,9*0,3,-6,3,0,-2,4,-2,10*0,-3,3,2*0,2,-2,2*0,&
                -1,1,6*0,3,-3,2*0,-2,2,5*0,1,-2,1,0,-2,4,-2,0,1,-2,1,9*0,-1,2,&
                -1,0,1,-2,1,10*0,1,-1,2*0,-1,1,6*0,-1,1,2*0,2,-2,2*0,-1,1/

        x(1:4)=y
        x(5:8)=y1*d1
        x(9:12)=y2*d2
        x(13:16)=y12*d1*d2
        x=matmul(wt,x)
        c=reshape(x,(/4,4/),order=(/2,1/))

end subroutine bcucof


!=================================================================================
        subroutine bcuint(y,y1,y2,y12,x1l,x1u,x2l,x2u,x1,x2,ansy,ansy1,ansy2)
!=================================================================================
        implicit none
        integer :: i
	real(kind=8), dimension(4), intent(in) :: y,y1,y2,y12
	real(kind=8), intent(in) :: x1l,x1u,x2l,x2u,x1,x2
	real(kind=8), intent(out) :: ansy,ansy1,ansy2
	real(kind=8) :: t,u
	real(kind=8), dimension(4,4) :: c

        call bcucof(y,y1,y2,y12,x1u-x1l,x2u-x2l,c)
        if (x1u == x1l .or. x2u == x2l) then
           print *, 'bcuint: problem with input values - boundary pair equal?'
        endif  
        t=(x1-x1l)/(x1u-x1l)
        u=(x2-x2l)/(x2u-x2l)
        ansy=0d0
        ansy2=0d0
        ansy1=0d0
        do i=4,1,-1
                ansy=t*ansy+((c(i,4)*u+c(i,3))*u+c(i,2))*u+c(i,1)
                ansy2=t*ansy2+(3d0*c(i,4)*u+2d0*c(i,3))*u+c(i,2)
                ansy1=u*ansy1+(3d0*c(4,i)*t+2d0*c(3,i))*t+c(2,i)
        end do
        ansy1=ansy1/(x1u-x1l)
        ansy2=ansy2/(x2u-x2l)

        end subroutine bcuint

  
!============================================================
 subroutine timeseries(n)
! Calculate and store global averages
!============================================================
   use share; use domain; use markeraux; use flow
   implicit none
   integer :: i,j,n
   real (kind=8) :: z,Nu_top,Nu_bot,Nu_top_min,Nu_top_max,Nu_bot_min,Nu_bot_max,    &
                    Vrms,Vrms_min,Vrms_max,Vsmax,Vsrms,Tav,Gav,etaav,Phiav,RaHav,   &
                    deplav,Phi_min,Phi_max,f1,f2,f3,f4,W,rho_r,drho_r,T_r,diff,eff, &
                    f_ppv,g,gmin,gmax,mixfrac,lid_depth
   real (kind=8), dimension(size(delz)) :: hflav,vpr
!===============================================================================================

!----------------------------------------------
! Top and bottom Nusselt numbers:
!----------------------------------------------
   Nu_top=0d0; Nu_bot=0d0
   Nu_top_min=1d50; Nu_top_max=-1d50
   Nu_bot_min=1d50; Nu_bot_max=-1d50
   do i=1,nx+1

     f1=(Tnew(i-1,nz+1)-Tnew(i-1,nz))/delz(nz+1)
     f2=(Tnew(i,nz+1)-Tnew(i,nz))/delz(nz+1)
     Nu_top=Nu_top-delx(i)/2d0*(f1+f2)
     if(f2.le.Nu_top_min) Nu_top_min=f2
     if(f2.ge.Nu_top_max) Nu_top_max=f2

     f1=(Tnew(i-1,1)-Tnew(i-1,0))/delz(1)
     f2=(Tnew(i,1)-Tnew(i,0))/delz(1)
     Nu_bot=Nu_bot-delx(i)/2d0*(f1+f2)
     if(f2.le.Nu_bot_min) Nu_bot_min=f2
     if(f2.ge.Nu_bot_max) Nu_bot_max=f2

   end do
   Nu_top=Nu_top/dimx
   Nu_bot=Nu_bot/dimx
   
!-----------------------------------------------------
! RMS velocity 'Vrms', max surface velocity 'Vsmax' 
! and rms surface velocity 'Vsrms'
!-----------------------------------------------------
   Vrms=0d0 
   Vrms_min=1d50; Vrms_max=-1d50
   do i=1,nx+1
     do j=1,nz+1

       f1=u(i-1,j-1)**2+v(i-1,j-1)**2
       f2=u(i,j-1)**2+v(i,j-1)**2
       f3=u(i-1,j)**2+v(i-1,j)**2
       f4=u(i,j)**2+v(i,j)**2
       Vrms=Vrms+delx(i)/2*delz(j)/2*(f1+f2+f3+f4)

       if(f4.le.Vrms_min) Vrms_min=f4
       if(f4.ge.Vrms_max) Vrms_max=f4
       
     end do
   end do
   Vrms=dsqrt(Vrms)/dimx/dimz
   Vrms_min=dsqrt(Vrms_min)/dimx/dimz
   Vrms_max=dsqrt(Vrms_max)/dimx/dimz

   Vsmax=maxval(dabs(u(:,nz+1)))
   
   Vsrms=0d0
   do i=1,nx+1
     Vsrms=Vsrms+delx(i)/2*(u(i-1,nz+1)**2 + u(i,nz+1)**2)
   end do
   Vsrms=dsqrt(Vsrms)/dimx

!----------------------------------
! Thickness of the stagnant lid
!----------------------------------
! Caluclate velocity profile
   vpr=0d0
   do j=1,nz+1
     do i=1,nx+1

       f1=u(i-1,j)**2+v(i-1,j)**2
       f2=u(i,j)**2+v(i,j)**2
       vpr(j)=vpr(j)+delx(i)/2d0*(f1+f2)

     end do
   end do
   vpr=dsqrt(vpr)/dimx
   
   lid_depth=0d0
! Stagnant-lid depth as the depth above which the average velocity is a fraction
! 'lidfrac' (typically 0.01) of the maximum of the velocity profile
   do j=1,nz+1

     if(vpr(j).le.lidfac*maxval(vpr)) then
       lid_depth=dimz-zgrid(j)
       exit
     endif

   end do

!-----------------------
! Average temperature 
!------------------------
   Tav=0d0
   do i=1,nx+1
     do j=1,nz+1
       
       if(therm.eq.1) then
         f1=Tnew(i-1,j-1)
         f2=Tnew(i,j-1)
         f3=Tnew(i-1,j)
         f4=Tnew(i,j)
         Tav=Tav+delx(i)/2*delz(j)/2*(f1+f2+f3+f4)/dimx/dimz
       endif

     end do
   end do

!-----------------------------
! Average viscosity
!-----------------------------
   etaav=0d0
   do i=1,nx+1
     do j=1,nz+1

       f1=eta(i-1,j-1)
       f2=eta(i,j-1)
       f3=eta(i-1,j)
       f4=eta(i,j)
       etaav=etaav+delx(i)/2*delz(j)/2*(f1+f2+f3+f4)/dimx/dimz

     end do
   end do

!-----------------------------
! Average internal heating
!-----------------------------
   RaHav=0d0
   do i=1,nx+1
     do j=1,nz+1

       f1=RaH(i-1,j-1)
       f2=RaH(i,j-1)
       f3=RaH(i-1,j)
       f4=RaH(i,j)
       RaHav=RaHav+delx(i)/2*delz(j)/2*(f1+f2+f3+f4)/dimx/dimz

     end do
   end do

!----------------------------------------
! Average composition and depletion
!----------------------------------------
   Gav=0d0
   deplav=0d0
   if(markers.eq.1) then
     do i=1,nx+1
       do j=1,nz+1
          
         if(comp_mk.eq.1) then 
           f1=Gam(i-1,j-1)
           f2=Gam(i,j-1)
           f3=Gam(i-1,j)
           f4=Gam(i,j)
           Gav=Gav+delx(i)/2*delz(j)/2*(f1+f2+f3+f4)/dimx/dimz
         endif

         if(depl_mk.eq.1) then 
           f1=depl(i-1,j-1)
           f2=depl(i,j-1)
           f3=depl(i-1,j)
           f4=depl(i,j)
           deplav=deplav+delx(i)/2*delz(j)/2*(f1+f2+f3+f4)/dimx/dimz
         endif

       end do
     end do
   endif

!----------------------------------------
! Average viscous dissipation 'Phiav'
!----------------------------------------
   Phiav=0d0
   Phi_min=1d50; Phi_max=-1d50 
   do i=1,nx+1
     do j=1,nz+1

       f1=Phi(i-1,j-1)
       f2=Phi(i,j-1)
       f3=Phi(i-1,j)
       f4=Phi(i,j)
       if(Di.ne.0d0) then
         Phiav=Phiav+Di/Ra*delx(i)/2*delz(j)/2*(f1+f2+f3+f4)/dimx/dimz
         if(f4.le.Phi_min) Phi_min=f4*Di/Ra
         if(f4.ge.Phi_max) Phi_max=f4*Di/Ra
       else
         Phiav=Phiav+1d0/Ra*delx(i)/2*delz(j)/2*(f1+f2+f3+f4)/dimx/dimz
         if(f4.le.Phi_min) Phi_min=f4/Ra
         if(f4.ge.Phi_max) Phi_max=f4/Ra
       endif

     end do
   end do

!---------------------------------------------------------
! Work done in adiabatic volume changes
! hflav: horizontally averaged convective heat flux
! cfr. eqs (53) and (54) of Jarvis & McKenzie (JFM, 1980)
!---------------------------------------------------------
   hflav=0d0
   do j=1,nz+1
     z=zgrid(j)-delz(j)/2

     if(Di.ne.0d0) then 
       call refstate(z,rho_r,drho_r,T_r)

       do i=1,nx+1
         hflav(j)=hflav(j)+delx(i)/2*rho_r*                 &
                  (Tnew(i-1,j)+Tnew(i,j))*(v(i-1,j)+v(i,j)) &
                  +delx(i)/2*rho_r*Tsurf*(v(i-1,j)+v(i,j))
       end do

     else

       do i=1,nx+1
         hflav(j)=hflav(j)+delx(i)/2*(Tnew(i-1,j)+Tnew(i,j))*(v(i-1,j)+v(i,j)) 
       end do
     
     endif
   end do

   W=0d0
   do j=1,nz+1
     W=W+delz(j)/2*(hflav(j-1)+hflav(j))
   end do
! c'e' un qualche errore che non trovo, dividere per 2 alla fine
! fa tornare i conti
   if(Di.ne.0d0) then
     W=W*Di/dimx/dimz/2d0
   else  
     W=W/dimx/dimz/2d0
   endif  

   diff=(Phiav-W)/max(Phiav,W)*100d0

!------------------------------------------------------
! Efficiency 'eff':
! average viscous dissipation over bottom heat flux
!-------------------------------------------------------
   eff=0d0
   do i=1,nx+1
     f1=(Tnew(i-1,1)-Tnew(i-1,0))/delz(1)
     f2=(Tnew(i,1)-Tnew(i,0))/delz(1)
     eff=eff+delx(i)/2*(f1+f2)
   end do
   eff=-Phiav/eff

! Mixed fraction (only meaningful for thermochemical convection)
! mixdepth is the depth below which the mixed fraction is to be computed
! (typically one has to look below the stagnant lid)
   gmin=10d0; gmax=-10d0
   do j=0,nz+1
     z=dimz-zgrid(j)

     if(z.ge.mixdepth) then
       g=0d0
       do i=1,nx+1
         g=g+delx(i)/2d0*(Gam(i,j)+Gam(i-1,j))/dimx
       end do
     endif

     if(g.gt.gmax) gmax=g
     if(g.lt.gmin) gmin=g

   end do
   mixfrac=dabs(gmin/gmax)

!-------------------------------------------
! PPV volume fraction
! f=(LD - int(G_pvppv dxdz) ) / LD
! LD is the domain volume and int is the
! integral of the ppv phase function
!-------------------------------------------
   f_ppv=0d0
!   if(pv_ppv.ne.0) then  

!   do i=1,nx+1
!     do j=1,nz+1

!       f1=phase_f(4,i-1,j-1)
!       f2=phase_f(4,i,j-1)
!       f3=phase_f(4,i-1,j)
!       f4=phase_f(4,i,j)
!       f_ppv=f_ppv+(f1+f2+f3+f4)*delx(i)/2*delz(j)/2
!
!     end do
!   end do
!   f_ppv=1d0-(dimx*dimz-f_ppv)/(dimx*dimz)

   open(1,file=trim(outdir)//'timeseries.dat',position='append')
   if(ti.eq.0) write(1,79) 'time', 'Nu', '<T>', 'Vrms', '<Vsurf>', '<eta>', '<Phi>', '<C>', 'Nu_min', 'Nu_max', &
                           'Vrms_min', 'Vrms_max', 'Phi_min', 'Phi_max', 'f_ppv', 'Mix_fraction', '<RaH>',      &
                           'Depletion', 'Lid-depth'
   write(1,80) ntime, Nu_top, Tav, Vrms, Vsrms, etaav, Phiav, Gav, Nu_top_min, Nu_top_max,      &
                 Vrms_min, Vrms_max, Phi_min, Phi_max, f_ppv, mixfrac, RaHav, deplav, lid_depth
   close(1)

   write(*,81) 'Time step=', n,'Time=',ntime,'<T>=',Tav,'<V>=',Vrms,'<Nu_top>=',Nu_top,'<Nu_bot>=',Nu_bot,'<eta>=',etaav,&
                '<C>=',Gav,'<F>=',deplav
!   write(*,81) 'Time step=', n,'Time=',ntime,'<T>=',Tav,'<V>=',Vrms,'<Nu>=',Nu_top,'<eta>=',etaav,'<C>=',Gav,'<F>=',deplav


79 format(20A15)
80 format(20E15.7)
81 format(A10,I8,A6,E12.5,A5,E12.5,A5,E12.5,A10,E12.5,A10,E12.5,A7,E12.5,A5,E12.5,A5,E12.5)


 end subroutine timeseries


!============================================================
 subroutine calc_profiles
! Calculate horizontally averaged profiles
!============================================================
   use share; use domain; use markeraux; use flow
   implicit none
   integer :: i,j
   real (kind=8) :: a1,a2
!============================================================

! Temperature
   T_z=0d0
   do j=0,nz+1
     do i=1,nx+1
       T_z(j)=T_z(j)+delx(i)/2d0*(Tnew(i,j)+Tnew(i-1,j))
     end do
   end do
   T_z=T_z/dimx
   
! Velocities
   v_z=0d0
   ux_z=0d0
   uz_z=0d0
   do j=0,nz+1
     do i=1,nx+1

     ! RMS
       a1=u(i-1,j)**2+v(i-1,j)**2
       a2=u(i,j)**2+v(i,j)**2
       v_z(j)=v_z(j)+delx(i)/2d0*(a1+a2)

     ! Vertical 
       uz_z(j)=uz_z(j)+delx(i)/2d0*(v(i,j)**2+v(i-1,j)**2)

     ! Horizontal
       ux_z(j)=ux_z(j)+delx(i)/2d0*(u(i,j)**2+u(i-1,j)**2)

     end do
   end do
   v_z=dsqrt(v_z/dimx)
   ux_z=dsqrt(ux_z/dimx)
   uz_z=dsqrt(uz_z/dimx)
   
! Viscosity
   eta_z=0d0
   do j=0,nz+1
     do i=1,nx+1
       eta_z(j)=eta_z(j)+delx(i)/2d0*(eta(i,j)+eta(i-1,j))
     end do
   end do
   eta_z=eta_z/dimx

! Viscous dissipation
   Phi_z=0d0
   do j=0,nz+1
     do i=1,nx+1
       Phi_z(j)=Phi_z(j)+delx(i)/2d0*(Phi(i,j)+Phi(i-1,j))
     end do
   end do
   if(Di.ne.0d0) then
      Phi_z=Di/Ra*Phi_z/dimx
   else   
      Phi_z=1d0/Ra*Phi_z/dimx
   endif
   
   if(markers.eq.1) then
!    Composition
     if(comp_mk.eq.1) then
       Gam_z=0d0
       do j=0,nz+1
         do i=1,nx+1
           Gam_z(j)=Gam_z(j)+delx(i)/2d0*(Gam(i,j)+Gam(i-1,j))
         end do
       end do
       Gam_z=Gam_z/dimx
     endif

! Depletion
     if(depl_mk.eq.1) then
       depl_z=0d0
       do j=0,nz+1
         do i=1,nx+1
           depl_z(j)=depl_z(j)+delx(i)/2d0*(depl(i,j)+depl(i-1,j))
         end do
       end do
       depl_z=depl_z/dimx
     endif
   endif

! Solidus
   Tsol_z=0d0
   do j=0,nz+1
     do i=1,nx+1
       Tsol_z(j)=Tsol_z(j)+delx(i)/2d0*(Tsol(i,j)+Tsol(i-1,j))
     end do
   end do
   Tsol_z=Tsol_z/dimx

! Heat sources
   RaH_z=0d0
   do j=0,nz+1
     do i=1,nx+1
       RaH_z(j)=RaH_z(j)+delx(i)/2d0*(RaH(i,j)+RaH(i-1,j))
     end do
   end do
   RaH_z=RaH_z/dimx

! Expansivity
   alpha_z=0d0
   do j=0,nz+1
     do i=1,nx+1
       alpha_z(j)=alpha_z(j)+delx(i)/2d0*(alpha(i,j)+alpha(i-1,j))
     end do
   end do
   alpha_z=alpha_z/dimx

! Conductivity
   kap_z=0d0
   do j=0,nz+1
     do i=1,nx+1
       kap_z(j)=kap_z(j)+delx(i)/2d0*(kap(i,j)+kap(i-1,j))
     end do
   end do
   kap_z=kap_z/dimx


 end subroutine calc_profiles


!============================================================
 subroutine write_output
! Write output files
!============================================================
   use share; use domain; use markeraux; use flow
   implicit none
   integer :: i,j,k
   real (kind=8) :: x,z,a
   real (kind=8), dimension(5) :: gfunc
   character (len=6) :: ci
   character (len=25) :: grd,profiles,mks
!============================================================

   timeref=ntime
   
   if(ti.lt.10) then
     write(ci,"('00000'I1)") ti
   elseif(ti.ge.10.and.ti.lt.100) then
     write(ci,"('0000'I2)") ti
   elseif(ti.ge.100.and.ti.lt.1000) then
     write(ci,"('000'I3)") ti
   elseif(ti.ge.1000.and.ti.lt.10000) then
     write(ci,"('00'I4)") ti
   elseif(ti.ge.10000.and.ti.lt.100000) then
     write(ci,"('0'I5)") ti
   else
     write(ci,"(I6)") ti
   endif

   grd="grd_"//trim(ci)//"_.dat"
   open(1,file=trim(outdir)//grd)
   
   open(2,file=trim(outdir)//'grd_restart.dat')
   
   profiles="profiles_"//trim(ci)//"_.dat"
   open(3,file=trim(outdir)//profiles)

   if(markers.gt.0.and.write_mks.eq.1) then
! Interpolate temperature and velocity on markers
     do i=1,Nmk
       Tmk(i)=0d0
       call grid2mk(xmk(i),zmk(i),Tnew,a)
       Tmk(i)=a

       umk(i)=0d0
       call grid2mk(xmk(i),zmk(i),u,a)
       umk(i)=a

       vmk(i)=0d0
       call grid2mk(xmk(i),zmk(i),v,a)
       vmk(i)=a
     end do

     mks="mks_"//trim(ci)//"_.dat"
     open(4,file=trim(outdir)//mks)
   endif

!---------------
! Write output 
!---------------
   gfunc=0d0
   write(1,83) 'Time step=',ti,'Time=',ntime,'x-points=',nx+2,'z-points=',nz+2,'aspect-ratio=',dimx
   write(2,83) 'Time step=',ti,'Time=',ntime,'x-points=',nx+2,'z-points=',nz+2,'aspect-ratio=',dimx
   do i=0,nx+1
     x=xgrid(i)
     do j=0,nz+1
       z=zgrid(j)
       
       do k=1,nphases
         gfunc(k)=phase_f(k,i,j)
       end do

       write(1,80) x,z,u(i,j),v(i,j),urms(i,j),Tnew(i,j),eta(i,j),Phi(i,j),Gam(i,j),RaH(i,j), &
                   alpha(i,j),kap(i,j),gfunc(1),gfunc(2),gfunc(3),gfunc(4),gfunc(5),          &
                   Tsol(i,j),depl(i,j),epsII(i,j),visc_mech(i,j)
       write(2,80) x,z,u(i,j),v(i,j),urms(i,j),Tnew(i,j),eta(i,j),Phi(i,j),Gam(i,j),RaH(i,j), &
                   alpha(i,j),kap(i,j),gfunc(1),gfunc(2),gfunc(3),gfunc(4),gfunc(5),          &
                   Tsol(i,j),depl(i,j),epsII(i,j),visc_mech(i,j)
     end do
   end do

   write(3,83) 'Time step=',ti,'Time=',ntime,'x-points=',nx+2,'z-points=',nz+2,'aspect-ratio=',dimx
   do j=0,nz+1
     z=zgrid(j)
     write(3,80) z,v_z(j),T_z(j),eta_z(j),Phi_z(j),Gam_z(j),ux_z(j),uz_z(j),RaH_z(j),alpha_z(j),kap_z(j), &
                 depl_z(j),Tsol0(j),Tsol_z(j),Tliq(j)
   end do
   close(1)
   close(2)
   close(3)
  
   if(markers.gt.0.and.write_mks.eq.1) then
     write(4,84) 'Time step=',ti,'Time=',ntime,'x-points=',nx+2,'z-points=',nz+2,'aspect-ratio=',dimx,'Markers=',Nmk
     do i=1,Nmk
       write(4,81) i,xmk(i),zmk(i),Tmk(i),umk(i),vmk(i)
     end do
   endif
   close(4)

!  Compute velocity gradient tensor for anisotropy calculations
!     if(markers.eq.1.and.write_mks.eq.1) then
!       write(3,83) 'Time step=',ti,'Time=',ntime
!       do i=1,Nmk
!         call grid2mk(xmk(i),zmk(i),dudx,a1)
!         call grid2mk(xmk(i),zmk(i),dvdz,a2)
!         call grid2mk(xmk(i),zmk(i),dvdx,a3)
!         call grid2mk(xmk(i),zmk(i),dudz,a4)	   
!         write(3,88) i,mktype(i),xmk(i),zmk(i),a1,a2,a3,a4
!       end do
!       
!     close(3)
!     endif
!88 format(2I8,2E12.4,10E15.5)

!===========================  Henry's code to story binary output  ===================================
!   if(mod(ti,stp_bin).eq.0.and.binout.eq.1) then 
! velocity storage in binary direct access files
!    WRITE(*,*) '==> Storing vx  in :vx_bin.dat (binary), rec #=',ir
!    WRITE(101,REC=ir) SNGL(ntime),SNGL(u(0:nx+1,0:nz+1))
!    WRITE(*,*) '==> Storing vz  in :vz_bin.dat (binary), rec #=',ir
!    WRITE(102,REC=ir) SNGL(ntime),SNGL(v(0:nx+1,0:nz+1))
!   endif       
!=================================================================================================    
80 format(20E20.8,I3)
81 format(I8,20E20.8)
83 format(A10,I7,2X,A5,E18.8,2X,A10,I6,A10,I6,A14,F5.2)
84 format(A10,I7,2X,A5,E18.8,2X,A10,I6,A10,I6,A14,F5.2,A9,I8)


 end subroutine write_output




!#######################################
 end module energy
!####################################### 

